% ----------------------------------------------------------------------------------------------- %
%                                                                                                 %
%                                   COUNTERFACTUAL POLICY SIMULATIONS                             %
%                                                                                                 %
% ----------------------------------------------------------------------------------------------- %


% We consider counterfactual analysis for the following policy:
% Students with $T_i=1$ are given priority over those with $T_i=0$, while
% within each type they are still ranked according to their indices. That
% is, given $T_i=T_{i^{\prime}}$, $i$ is ranked higher by all schools than
% $i^{\prime}$ if and only if $i<i^{\prime}$. One may consider this as an
% affirmative action policy if $T_i=1$ indicates that $i$ is in a
% disadvantaged group. The algorithm, the random serial dictatorship, is
% kept the same.


% We carry out the following counterfactual analyses.
%
% (1) True Preferences with mistakes: We use the true coefficients in utility
%   functions to simulate counterfactual outcomes. Students submit truthful 
%   12-school rank-order lists. This is taken as a benchmark.
%
% (2) Submitted ROLs: One assumes that the submitted ROLs are true
% ordinal preferences and that students submit the same ROLs given the
% counterfactual policy.
%
% (3) Truth-telling Estimates: One assumes that the submitted ROLs are true
% ordinal preferences, and therefore student preferences can be estimated
% from the data with truth-telling as the identifying condition. With the
% counterfactual policy, we simulate student preferences based on the
% estimates and let students submit truthful 12-school rank-order lists.
%
% (4) Stability Estimates: We estimate student preferences from the data
% with stability as the identifying condition. With the counterfactual
% policy, we simulate student preferences based on the estimates and let
% students submit truthful 12-school rank-order lists.


% ----------------------------------------------------------------------------------------------- %
%%                                         PREPARE DATA                                            %
% ----------------------------------------------------------------------------------------------- %

tic
cd(folder_data)
load sim_data_counterfactual % simulated school choice data; no estimates

% counterfactual: new priority ranking
% MC.Priority = (1:I)'/I; % same priority at each school.
% Student binary characteristics
% rng(101);
% MC.Stu_bin = (rand(I,M)<1/3); % minorities
% MC.Stu_char = MC.Stu_bin;

CF.Priority = repmat(MC.Priority, [1,M]) + MC.Stu_char;
for mm = 1:M
    [~,ind] = sort(CF.Priority(:,mm),'ascend');
    
    CF.Priority_rank(ind,mm) = (1:I)/I;
end
% ------------------------- %
% 3 SETS OF POINT-ESTIMATES %
% ------------------------- %

% estimates when there are mistakes; the 1st is the TT case (no skips)
load est_Mis.mat LL_TT_est LL_ST_est
% load est_R1 LL_R1_est
% LL_Ro_est = LL_R1_est;
% ----------------------------------------------------------------------------------------------- %
%                                                                                                 %
%%                                  PART I : Simulate w/ True Pref (no mistake) or Submitted ROLs %
%                                                                                                 %
% ----------------------------------------------------------------------------------------------- %
% Results collected in CFROL

% ------------------------------------------------------ %
% Run Student-Proposing DA algorithm                     %
% ------------------------------------------------------ %

CF_Assigned=NaN(M,I,P);
CF_Matched=NaN(M,J,I,P);
CF_Cutoffs=NaN(J,M,P);

% Loop over MC samples
for mm = 1:M
    Priority_mm = CF.Priority_rank(:,mm);
    for pp = 1:P
        % Assigned: each student's assigned school
        % Matched: each school's assigned students
        [Assigned, Matched] = func_DA('stu', squeeze(Mis.Stu_rk_pref(:,:,mm,pp))', repmat(Priority_mm',J,1), Capacities);
        
        % Find school cutoffs
        Cutoffs=NaN(J,1);
        for jj = 1:J % loop over schools
            Cutoffs(jj) = min(Priority_mm(Matched(jj,:) == 1));
            % If school capacity is not filled, then cutoff is 0
            if sum(Matched(jj,:)) < Capacities(jj)
                Cutoffs(jj)=0;
            end
            % If the school's cutoff is the student with smallest priority, then cutoff is zero
            if Cutoffs(jj) == min(Priority_mm)
                Cutoffs(jj)=0;
            end
        end
        CF_Assigned(mm,:,pp) = Assigned;
        CF_Matched(mm,:,:,pp) = Matched;
        CF_Cutoffs(:,mm,pp) = Cutoffs;
    end
end

CFROL.Assigned = CF_Assigned;
CFROL.Matched = CF_Matched;
CFROL.Cutoffs = CF_Cutoffs;

clear CF_Assigned CF_Matched CF_Cutoffs

% ----------------------------------------------------------------------------------------------- %
%                                                                                                 %
%%                       PART II : Simulate with WTT Estimates                                     %
%                                                                                                 %
% ----------------------------------------------------------------------------------------------- %
% Results collected in CFWTT

% ------------------------------------------------------ %
% GENERATE PREFERENCES WITH TT ESTIMATES                 %
% ------------------------------------------------------ %

CFWTT.U_ij = NaN(I,J,M,P);
CFWTT.ROL = NaN(I,J,M,P);

for pp = 1:P
    for mm = 1:M
        coeff_est = LL_TT_est(:,(pp-1)*M+mm);
        CFWTT.U_ij(:,:,mm,pp) = repmat([1:J]*coeff_est(1),[I,1]) ... % repmat([0;coeff_est(1:J-1)]',[I,1]) ...
            + coeff_est(end-1) .* MC.distance_school(:,:,mm) ...
            + coeff_est(end-2) .* School_char(:,:,mm) .* repmat(MC.Stu_char(:,mm),[1,J]) ...
            + coeff_est(end) .* School_small(:,:,mm) ...
            + MC.E_ij(:,:,mm) + 1000; % make it positive
%         CFWTT.U_ij(:,:,mm,pp) = repmat([1:J]*coeff_est(1),[I,1]) ...
%             - MC.distance_school(:,:,mm) ...
%             + coeff_est(end-1) .* School_char(:,:,mm) .* repmat(MC.Stu_char(:,mm),[1,J]) ...
%             + MC.E_ij(:,:,mm).*exp(coeff_est(end)) ;
    end
    [~, CFWTT.ROL(:,:,:,pp)] = sort(CFWTT.U_ij(:,:,:,pp),2,'descend');
end


% ------------------------------------------------------ %
% Run Student-Proposing DA algorithm                     %
% ------------------------------------------------------ %


CF_Assigned=NaN(M,I,P);
CF_Matched=NaN(M,J,I,P);
CF_Cutoffs=NaN(J,M,P);

% Loop over MC samples
for mm = 1:M
        Priority_mm = CF.Priority_rank(:,mm);
    for pp = 1:P
        % Assigned: each student's assigned school
        % Matched: each school's assigned students
        [Assigned, Matched] = func_DA('stu', squeeze(CFWTT.U_ij(:,:,mm,pp))', repmat(Priority_mm',J,1), Capacities);
        
        % Find school cutoffs
        Cutoffs=NaN(J,1);
        for jj = 1:J % loop over schools
            Cutoffs(jj) = min(Priority_mm(Matched(jj,:) == 1));
            % If school capacity is not filled, then cutoff is 0
            if sum(Matched(jj,:)) < Capacities(jj)
                Cutoffs(jj)=0;
            end
            % If the school's cutoff is the student with smallest priority, then cutoff is zero
            if Cutoffs(jj) == min(Priority_mm)
                Cutoffs(jj)=0;
            end
        end
        CF_Assigned(mm,:,pp) = Assigned;
        CF_Matched(mm,:,:,pp) = Matched;
        CF_Cutoffs(:,mm,pp) = Cutoffs;
    end
end

CFWTT.Assigned = CF_Assigned;
CFWTT.Matched = CF_Matched;
CFWTT.Cutoffs = CF_Cutoffs;

clear CF_Assigned CF_Matched CF_Cutoffs


% ----------------------------------------------------------------------------------------------- %
%                                                                                                 %
%%                       PART III : Simulate with STability Estimates                             %
%                                                                                                 %
% ----------------------------------------------------------------------------------------------- %
% Results collected in CFST

% ------------------------------------------------------ %
% GENERATE PREFERENCES WITH TT ESTIMATES                 %
% ------------------------------------------------------ %

CFST.U_ij = NaN(I,J,M,P);
CFST.ROL = NaN(I,J,M,P);

for pp = 1:P
    for mm = 1:M
        coeff_est = LL_ST_est(:,(pp-1)*M+mm);
        CFST.U_ij(:,:,mm,pp) = repmat([1:J]*coeff_est(1),[I,1]) ... % repmat([0;coeff_est(1:J-1)]',[I,1]) ...
            + coeff_est(end-1) .* MC.distance_school(:,:,mm) ...
            + coeff_est(end-2) .* School_char(:,:,mm) .* repmat(MC.Stu_char(:,mm),[1,J]) ...
            + coeff_est(end) .* School_small(:,:,mm) ...
            + MC.E_ij(:,:,mm) + 1000; % make it positive
%         CFST.U_ij(:,:,mm,pp) = repmat([1:J]*coeff_est(1),[I,1]) ...
%             - MC.distance_school(:,:,mm) ...
%             + coeff_est(end-1) .* School_char(:,:,mm) .* repmat(MC.Stu_char(:,mm),[1,J]) ...
%             + MC.E_ij(:,:,mm).*exp(coeff_est(end)) ;        
    end
    [~, CFST.ROL(:,:,:,pp)] = sort(CFST.U_ij(:,:,:,pp),2,'descend');
end


% ------------------------------------------------------ %
% Run Student-Proposing DA algorithm                     %
% ------------------------------------------------------ %


CF_Assigned=NaN(M,I,P);
CF_Matched=NaN(M,J,I,P);
CF_Cutoffs=NaN(J,M,P);

% Loop over MC samples
for mm = 1:M
        Priority_mm = CF.Priority_rank(:,mm);    
    for pp = 1:P
        % Assigned: each student's assigned school
        % Matched: each school's assigned students
        [Assigned, Matched] = func_DA('stu', squeeze(CFST.U_ij(:,:,mm,pp))', repmat(Priority_mm',J,1), Capacities);
        
        % Find school cutoffs
        Cutoffs=NaN(J,1);
        for jj = 1:J % loop over schools
            Cutoffs(jj) = min(Priority_mm(Matched(jj,:) == 1));
            % If school capacity is not filled, then cutoff is 0
            if sum(Matched(jj,:)) < Capacities(jj)
                Cutoffs(jj)=0;
            end
            % If the school's cutoff is the student with smallest priority, then cutoff is zero
            if Cutoffs(jj) == min(Priority_mm)
                Cutoffs(jj)=0;
            end
        end
        CF_Assigned(mm,:,pp) = Assigned;
        CF_Matched(mm,:,:,pp) = Matched;
        CF_Cutoffs(:,mm,pp) = Cutoffs;
    end
end

CFST.Assigned = CF_Assigned;
CFST.Matched = CF_Matched;
CFST.Cutoffs = CF_Cutoffs;

clear CF_Assigned CF_Matched CF_Cutoffs


% ----------------------------------------------------------------------------------------------- %
%                                                                                                 %
%%                                  PART IV: Simulate w/ True Pref (w/ mistake)                   %
%                                                                                                 %
% ----------------------------------------------------------------------------------------------- %
% mistake_pp: mistake percentages
% use the same preference profiles to simulate cutoff distribution

% prob(making a mistake is positively correlated with Stu_char)

% prob(making a mistake is positively correlated with Stu_char)
rng(106);
rand1 = rand(I,M);

Mis.Mistake = (rand1 > 0.2);
Mis.Mistake = repmat(Mis.Mistake,[1,1,P]);

% ---------------------------------- %
% RUN DA TO DETERMINE SCHOOL CUTOFFS %
% ---------------------------------- %

% priorities
CDCF.Priority = repmat(CD.Priority, [1,M_CD]) + CD.Stu_char;

for mm = 1:M_CD
    [~,ind] = sort(CDCF.Priority(:,mm),'ascend');
    CDCF.Priority_rank(ind,mm) = (1:I)/I;
end

CDCF.Cutoffs = NaN(J,M_CD);
CDCF.Assigned = NaN(M_CD,I);
CDCF.Matched =  NaN(M_CD,J,I);

% Run Student-Proposing DA algorithm

CDCF_Assigned=NaN(M_CD,I);
CDCF_Matched=NaN(M_CD,J,I);
CDCF_Cutoffs=NaN(J,M_CD);


% Loop over CD (cutoff distribution) samples
parfor mm = 1:M_CD
    Priority_mm = CDCF.Priority_rank(:,mm);
    % Assigned: each student's assigned school
    % Matched: each school's assigned students
    [Assigned, Matched] = func_DA('stu', CD.Stu_rk_pref(:,:,mm)', repmat(Priority_mm',J,1), Capacities);
    
    % Find school cutoffs
    Cutoffs=NaN(J,1);
    for jj = 1:J % loop over schools
        Cutoffs(jj) = min(Priority_mm(Matched(jj,:) == 1) );
        % If school capacity is not filled, then cutoff is 0
        if sum(Matched(jj,:)) < Capacities(jj)
            Cutoffs(jj)=0;
        end
        % If the school's cutoff is the student with smallest priority, then cutoff is zero
        if Cutoffs(jj) == 1/I
            Cutoffs(jj)=0;
        end
    end
    CDCF_Assigned(mm,:) = Assigned;
    CDCF_Matched(mm,:,:) = Matched;
    CDCF_Cutoffs(:,mm) = Cutoffs;
end

CDCF.Assigned = CDCF_Assigned;
CDCF.Matched = CDCF_Matched;
CDCF.Cutoffs = CDCF_Cutoffs;

clear CD_Assigned CD_Matched CD_Cutoffs;

% admission probability when truth-telling

% Probability of being accepted by each school given the true preference
% order
% ROL: (IxJ) matrix where (i,j) = identifier of school ranked j-th in student i's preferences (from top = col 1 to bottom = col J)

TTCF.PA = NaN(I,J,M);
TTCF.PA_reorder = TTCF.PA;
TTCF.PA_ro = TTCF.PA;
% Loop over MC samples
for ii = 1:I
    for mm=1:M
            % prob being accepted by one's first choice
            TTCF.PA(ii,1,mm) = sum(CF.Priority_rank(ii,mm) >= CDCF.Cutoffs(TT.ROL(ii,1,mm),:))/M_CD;
            
            TTCF.PA(ii,2,mm) = sum((CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,1,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) >= CDCF.Cutoffs(TT.ROL(ii,2,mm),:)))/M_CD;
            
            TTCF.PA(ii,3,mm) = sum((CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,1,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,2,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) >= CDCF.Cutoffs(TT.ROL(ii,3,mm),:)))/M_CD;
            
            TTCF.PA(ii,4,mm) = sum((CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,1,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,2,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,3,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) >= CDCF.Cutoffs(TT.ROL(ii,4,mm),:)))/M_CD;
            
            TTCF.PA(ii,5,mm) = sum((CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,1,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,2,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,3,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,4,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) >= CDCF.Cutoffs(TT.ROL(ii,5,mm),:)))/M_CD;
            
            TTCF.PA(ii,6,mm) = sum((CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,1,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,2,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,3,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,4,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,5,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) >= CDCF.Cutoffs(TT.ROL(ii,6,mm),:)))/M_CD;
            
            TTCF.PA(ii,7,mm) = sum((CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,1,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,2,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,3,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,4,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,5,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,6,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) >= CDCF.Cutoffs(TT.ROL(ii,7,mm),:)))/M_CD;
            
            TTCF.PA(ii,8,mm) = sum((CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,1,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,2,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,3,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,4,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,5,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,6,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,7,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) >= CDCF.Cutoffs(TT.ROL(ii,8,mm),:)))/M_CD;
            
            TTCF.PA(ii,9,mm) = sum((CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,1,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,2,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,3,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,4,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,5,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,6,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,7,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,8,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) >= CDCF.Cutoffs(TT.ROL(ii,9,mm),:)))/M_CD;
            
            TTCF.PA(ii,10,mm) = sum((CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,1,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,2,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,3,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,4,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,5,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,6,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,7,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,8,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,9,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) >= CDCF.Cutoffs(TT.ROL(ii,10,mm),:)))/M_CD;
            
            TTCF.PA(ii,11,mm) = sum((CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,1,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,2,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,3,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,4,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,5,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,6,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,7,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,8,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,9,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,10,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) >= CDCF.Cutoffs(TT.ROL(ii,11,mm),:)))/M_CD;
            
            TTCF.PA(ii,12,mm) = sum((CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,1,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,2,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,3,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,4,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,5,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,6,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,7,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,8,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,9,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,10,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) < CDCF.Cutoffs(TT.ROL(ii,11,mm),:)) ...
                .*(CF.Priority_rank(ii,mm) >= CDCF.Cutoffs(TT.ROL(ii,12,mm),:)))/M_CD;
            
            TTCF.PA_ro(ii,TT.ROL(ii,:,mm),mm) = TTCF.PA(ii,:,mm);
    end
end

TTCF.PA_reorder = TTCF.PA_ro;

% rng(50000)
% TT_ran = rand(I,1); % to be used in assigning the mistakes

% TTCF.PA_reorder((TTCF.PA_ro == 0) .* repmat((TT_ran <= 1/3),[1,J,M]) == 1) = -2.5;
% TTCF.PA_reorder((TTCF.PA_ro == 0) .* repmat((TT_ran > 1/3),[1,J,M]) .* repmat((TT_ran <= 2/3),[1,J,M]) == 1) = -1.5;
% TTCF.PA_reorder((TTCF.PA_ro == 0) .* repmat((TT_ran > 2/3),[1,J,M]) .* repmat((TT_ran <= 3/3),[1,J,M]) == 1) = -0.5;

% mistake_pp = [-10,-1.5,-0.5,0.2]'; % refer to PA_dis; new version with 4 DGPs: TT, IRR1-2, REL
rng(50000)
TT_ran = rand(I,J,M); % to be used in assigning the mistakes
TTCF.PA_reorder(((TTCF.PA_ro == 0) .* (TT_ran <= 1/4)) == 1) = -2;
TTCF.PA_reorder(((TTCF.PA_ro == 0) .* (TT_ran >  1/4).*(TT_ran <= 5/6) == 1)) = -2;
%TT.PA_reorder(((TT.PA_ro == 0) .* (TT_ran >  2/4).*(TT_ran <= 3/4) == 1)) = -2;
TTCF.PA_reorder(((TTCF.PA_ro == 0) .* (TT_ran >  5/6) == 1)) = 2;

clear TT_ran

% for those who have all admission probabilities = 0, keep a random number
% of them
rng(50001)
TT_ind = 1 - (0-floor(12*rand(I,M,12))).*(rand(I,M,12)>0.8);

for mm = 1:M
    for ii=1:I
        if sum(TTCF.PA_ro(ii,:,mm)) == 0
            %TT.PA_reorder(ii,TT.ROL(ii,1:2,mm),mm) = 2;
            TTCF.PA_reorder(ii,TT.ROL(ii,TT_ind(ii,mm,:),mm),mm) = 2;
        elseif max(TTCF.PA_ro(ii,:,mm)) > 0 && max(TTCF.PA_ro(ii,:,mm)) <= mistake_pp(P) % max admission prob to be skipped
            %TT.PA_reorder(ii,TT.ROL(ii,1:2,mm),mm) = 2;
            TTCF.PA_reorder(ii,TT.ROL(ii,TT_ind(ii,mm,:),mm),mm) = 2;
        end
    end
end

clear TT_ind


CFMis.ROL = NaN(I,J,M,P);
CFMis.Stu_rk_pref = NaN(I,J,M,P);
for pp = 1:P
    for jj = 1:J
        CFMis.Stu_rk_pref(:,jj,:,pp) = ...
            TT.Stu_rk_pref(:,jj,:).*(1 - ...
             (reshape(Mis.Mistake(:,:,pp),[I,1,M]) == 1) ...
            .* (TTCF.PA_reorder(:,jj,:) <= mistake_pp(pp)));
    end
    [u_val, CFMis.ROL(:,:,:,pp)] = sort(Mis.Stu_rk_pref(:,:,:,pp),2,'descend');
    CFMis.ROL(:,:,:,pp) = CFMis.ROL(:,:,:,pp) .* (u_val>0);
end


% ------------------------------------------------------ %
% Run Student-Proposing DA algorithm                     %
% ------------------------------------------------------ %

CFMis_Assigned=NaN(M,I,P);
CFMis_Matched=NaN(M,J,I,P);
CFMis_Cutoffs=NaN(J,M,P);

% Loop over MC samples
for mm = 1:M
    Priority_mm = CF.Priority_rank(:,mm);
    for pp = 1:P
        % Assigned: each student's assigned school
        % Matched: each school's assigned students
        [Assigned, Matched] = func_DA('stu', squeeze(CFMis.Stu_rk_pref(:,:,mm,pp))', repmat(Priority_mm',J,1), Capacities);
        
        % Find school cutoffs
        Cutoffs=NaN(J,1);
        for jj = 1:J % loop over schools
            Cutoffs(jj) = min(Priority_mm(Matched(jj,:) == 1));
            % If school capacity is not filled, then cutoff is 0
            if sum(Matched(jj,:)) < Capacities(jj)
                Cutoffs(jj)=0;
            end
            % If the school's cutoff is the student with smallest priority, then cutoff is zero
            if Cutoffs(jj) == min(Priority_mm)
                Cutoffs(jj)=0;
            end
        end
        CFMis_Assigned(mm,:,pp) = Assigned;
        CFMis_Matched(mm,:,:,pp) = Matched;
        CFMis_Cutoffs(:,mm,pp) = Cutoffs;
    end
end

CFMis.Assigned = CFMis_Assigned;
CFMis.Matched = CFMis_Matched;
CFMis.Cutoffs = CFMis_Cutoffs;

clear CFMis_Assigned CFMis_Matched CFMis_Cutoffs


% ----------------------------------------------------------------------------------------------- %
%                                                                                                 %
%%                                   PART V : CUTOFFS                                            %
%                                                                                                 %
% ----------------------------------------------------------------------------------------------- %

CF.mean_cutoff = NaN(J,4,P);

CF.mean_cutoff(:,1,:) = mean(CFROL.Cutoffs,2);
CF.mean_cutoff(:,2,:) = mean(CFWTT.Cutoffs,2);
CF.mean_cutoff(:,3,:) = mean(CFST.Cutoffs,2);
CF.mean_cutoff(:,4,:) = mean(CFMis.Cutoffs,2);

% CF.diff_cutoff = CF.mean_cutoff - repmat(CF.mean_cutoff(:,1,1),[1,4,P]);
CF.diff_cutoff = CF.mean_cutoff - repmat(CF.mean_cutoff(:,4,:),[1,4,1]);

% plot mean differences between true and predicted cutoffs.
% 8 graphs; in each graph, one ture, one ROL, one TT est, one ST est
formatSpec = 'fig_3_CF_cutoffs_%d';
figSpec = 'subfig (%d)';

cd(folder_figures)

for pp = 1:P
    figname = sprintf(figSpec, pp);
    figure()
    x = 1:J;
    p1 = plot(x,CF.diff_cutoff(:,1,pp),'-.g','LineWidth',4);
    hold on
    p2 = plot(x,CF.diff_cutoff(:,2,pp),'-r','LineWidth',4);
    hold on
    p3 = plot(x,CF.diff_cutoff(:,3,pp),':b','LineWidth',4);
    hold on    
    line([1 12], [0 0],'Color','black','LineWidth',0.5)

    %title(figname)
    ax = gca; % current axes
    ax.FontSize = 25;
    ax.TickDir = 'out';
    ax.TickLength = [0.01 0.01];
    ax.YLim = [-0.2 0.1];
    ax.YTick = ([-0.2 -0.1 0 0.1]);
    ax.XLim = [1 12];
    ax.XTick = ([2 4 6 8 10 12]);
    xlabel('College Index');
    if pp == 1
        l = legend([p1,p2,p3],{'Estimates w/ Submitted ROLs','WTT Estimates', 'Stability Estimates'},'Box','off','Location','southwest');
        set(l,'FontSize',20);
        ylabel('( Predicted Cutoff ) - ( True Cutoff )','FontSize',22);
    end    
    %     ax.XTick = ([1 2 3 4 5 6 7 8 9 10 11 12]);
    %     xticklabels({'-3\pi','-2\pi','-\pi','0','\pi','2\pi','3\pi'})
    filename = sprintf(formatSpec, pp);
    print(filename,'-djpeg')
    print(filename,'-depsc2')
end
close all

% ----------------------------------------------------------------------------------------------- %
%                                                                                                 %
%%                                   PART VII: MISPREDICTED MATCH                                 %
%                                                                                                 %
% ----------------------------------------------------------------------------------------------- %


% --------------- %
% Different Match %
% --------------- %

CF.Diff_Assign = NaN(3,P);
for pp = 1:P
    CF.Diff_Assign(1,pp) = mean(mean(CFROL.Assigned(:,:,pp) ~= CFMis.Assigned(:,:,pp)));
    CF.Diff_Assign(2,pp) = mean(mean(CFWTT.Assigned(:,:,pp) ~= CFMis.Assigned(:,:,pp)));
    CF.Diff_Assign(3,pp) = mean(mean(CFST.Assigned(:,:,pp) ~= CFMis.Assigned(:,:,pp)));
end

CF.Diff_Assign_AA = NaN(3,P); % for the minority group
for pp = 1:P
    CF.Diff_Assign_AA(1,pp) = sum(sum((CFROL.Assigned(:,:,pp) ~= CFMis.Assigned(:,:,pp))' .* MC.Stu_char)) / sum(sum(MC.Stu_char));
    CF.Diff_Assign_AA(2,pp) = sum(sum((CFWTT.Assigned(:,:,pp) ~= CFMis.Assigned(:,:,pp))' .* MC.Stu_char)) / sum(sum(MC.Stu_char));
    CF.Diff_Assign_AA(3,pp) = sum(sum((CFST.Assigned(:,:,pp) ~= CFMis.Assigned(:,:,pp))' .* MC.Stu_char)) / sum(sum(MC.Stu_char));
end

CF.Diff_Assign_nAA = NaN(3,P); % for the non-minority group
for pp = 1:P
    CF.Diff_Assign_nAA(1,pp) = sum(sum((CFROL.Assigned(:,:,pp) ~= CFMis.Assigned(:,:,pp))' .* (1-MC.Stu_char))) / sum(sum(1-MC.Stu_char));
    CF.Diff_Assign_nAA(2,pp) = sum(sum((CFWTT.Assigned(:,:,pp) ~= CFMis.Assigned(:,:,pp))' .* (1-MC.Stu_char))) / sum(sum(1-MC.Stu_char));
    CF.Diff_Assign_nAA(3,pp) = sum(sum((CFST.Assigned(:,:,pp) ~= CFMis.Assigned(:,:,pp))' .* (1-MC.Stu_char))) / sum(sum(1-MC.Stu_char));
end


% ------ %
%   AA   %
% ------ %

figure()
% area([2;3;4],[1*ones(3,1)],'FaceColor',[0.9  0.9  0.9],'LineStyle','none');
hold on
x = 1:P;
p1 = plot(x,CF.Diff_Assign_AA(1,:),'-.g','LineWidth',2);
hold on
p2 = plot(x,CF.Diff_Assign_AA(2,:),'-r','LineWidth',2);
hold on
p3 = plot(x,CF.Diff_Assign_AA(3,:),':b','LineWidth',2);
hold on
line([1 4], [CF.Diff_Assign_AA(1,1) CF.Diff_Assign_AA(1,1)],'Color','black','LineWidth',0.5)
legend([p1,p2,p3],{'Estimates w/ Submitted ROLs','WTT Estimates', 'Stability Estimates'},'Box','off','Location','northwest')

%title(figname)
ax = gca; % current axes
ax.FontSize = 14;
ax.TickDir = 'out';
ax.TickLength = [0.01 0.01];
ax.YLim = [0,0.5];
% ax.XLim = [1 4];
% ax.XTick = ([1 2 3 4]);
% ax.XTickLabel = ({'TT','IRR1','IRR2','REL'});
ax.XLim = [1 3];
ax.XTick = ([1 2 3]);
ax.XTickLabel = ({'TT','PIM','PRM'});
ylabel('Fraction w/ Mis-predicted Match');
xlabel('Data Generating Processes');
print('fig_4_CF_assign_AA','-djpeg')
print('fig_4_CF_assign_AA','-depsc2')

% ------ %
% non-AA %
% ------ %

figure()
% area([2;3;4],[1*ones(3,1)],'FaceColor',[0.9  0.9  0.9],'LineStyle','none');
hold on
x = 1:P;
p1 = plot(x,CF.Diff_Assign_nAA(1,:),'-.g','LineWidth',2);
hold on
p2 = plot(x,CF.Diff_Assign_nAA(2,:),'-r','LineWidth',2);
hold on
p3 = plot(x,CF.Diff_Assign_nAA(3,:),':b','LineWidth',2);
hold on
legend([p1,p2,p3],{'Estimates w/ Submitted ROLs','WTT Estimates', 'Stability Estimates'},'Box','off','Location','northwest')

%title(figname)
ax = gca; % current axes
ax.FontSize = 14;
ax.TickDir = 'out';
ax.TickLength = [0.01 0.01];
ax.YLim = [0,0.5];
% ax.XLim = [1 4];
% ax.XTick = ([1 2 3 4]);
% ax.XTickLabel = ({'TT','IRR1','IRR2','REL'});
ax.XLim = [1 3];
ax.XTick = ([1 2 3]);
ax.XTickLabel = ({'TT','PIM','PRM'});
ylabel('Fraction w/ Mis-predicted Match');
xlabel('Data Generating Processes');
print('fig_E2_CF_assign_nAA','-djpeg')
print('fig_E2_CF_assign_nAA','-depsc2')

close all


% --------------- %
% Welfare change  %
% --------------- %
CF.Wel = NaN(5,P,2); % better off, worse off
clear U_ij
U_ij = MC.U_ij;
U_ij(:,13,:) = -100; % arbitrary negative number

Mis.Assigned(Mis.Assigned==0) = 13; % assignment under the old policy
CFROL.Assigned(CFROL.Assigned == 0) = 13;
CFWTT.Assigned(CFWTT.Assigned == 0) = 13;
CFST.Assigned(CFST.Assigned == 0) = 13;
CFMis.Assigned(CFMis.Assigned == 0) = 13;

CFWTT.U_ij(:,13,:,:) = -100;
CFST.U_ij(:,13,:,:) = -100;
CFRo.U_ij(:,13,:,:) = -100;

CFMis.U_ij = repmat(U_ij,[1,1,1,P]);

CFROL.U_ij = Mis.Stu_rk_pref;
CFROL.U_ij(:,13,:,:) = -100;

for pp = 1:P
    % worse off
    wel = zeros(I,M,4);
    for mm = 1:M
        for ii = 1:I
            wel(ii,mm,1) = (CFROL.U_ij(ii,squeeze(Mis.Assigned(mm,ii,1)),mm,pp) > CFROL.U_ij(ii,squeeze(CFROL.Assigned(mm,ii,pp)),mm,pp));
            wel(ii,mm,2) = (CFWTT.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) > CFWTT.U_ij(ii,squeeze(CFWTT.Assigned(mm,ii,pp)),mm,pp));
            wel(ii,mm,3) = (CFST.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) > CFST.U_ij(ii,squeeze(CFST.Assigned(mm,ii,pp)),mm,pp));
            wel(ii,mm,4) = (CFMis.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) > CFMis.U_ij(ii,squeeze(CFMis.Assigned(mm,ii,pp)),mm,pp));
        end
    end
    CF.Wel(1,pp,1) = mean(mean(wel(:,:,1)));
    CF.Wel(2,pp,1) = mean(mean(wel(:,:,2)));
    CF.Wel(3,pp,1) = mean(mean(wel(:,:,3)));
    CF.Wel(4,pp,1) = mean(mean(wel(:,:,4)));
    % better off
    wel = zeros(I,M,4);
    for mm = 1:M
        for ii = 1:I
            wel(ii,mm,1) = (CFROL.U_ij(ii,squeeze(Mis.Assigned(mm,ii,1)),mm,pp) < CFROL.U_ij(ii,squeeze(CFROL.Assigned(mm,ii,pp)),mm,pp));
            wel(ii,mm,2) = (CFWTT.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) < CFWTT.U_ij(ii,squeeze(CFWTT.Assigned(mm,ii,pp)),mm,pp));
            wel(ii,mm,3) = (CFST.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) < CFST.U_ij(ii,squeeze(CFST.Assigned(mm,ii,pp)),mm,pp));
            wel(ii,mm,4) = (CFMis.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) < CFMis.U_ij(ii,squeeze(CFMis.Assigned(mm,ii,pp)),mm,pp));
        end
    end
    CF.Wel(1,pp,2) = mean(mean(wel(:,:,1)));
    CF.Wel(2,pp,2) = mean(mean(wel(:,:,2)));
    CF.Wel(3,pp,2) = mean(mean(wel(:,:,3)));
    CF.Wel(4,pp,2) = mean(mean(wel(:,:,4)));
end

CF.Wel_AA = NaN(4,P,3); % for the minority group
CF.Wel_AA_sd = NaN(4,P,3); % standard deviation

for pp = 1:P
    % worse off
    wel = zeros(I,M,4);
    for mm = 1:M
        for ii = 1:I
            wel(ii,mm,1) = (CFROL.U_ij(ii,squeeze(Mis.Assigned(mm,ii,1)),mm,pp) > CFROL.U_ij(ii,squeeze(CFROL.Assigned(mm,ii,pp)),mm,pp)) .* MC.Stu_char(ii,mm);
            wel(ii,mm,2) = (CFWTT.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) > CFWTT.U_ij(ii,squeeze(CFWTT.Assigned(mm,ii,pp)),mm,pp)) .* MC.Stu_char(ii,mm);
            wel(ii,mm,3) = (CFST.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) > CFST.U_ij(ii,squeeze(CFST.Assigned(mm,ii,pp)),mm,pp)) .* MC.Stu_char(ii,mm);
            wel(ii,mm,4) = (CFMis.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) > CFMis.U_ij(ii,squeeze(CFMis.Assigned(mm,ii,pp)),mm,pp)) .* MC.Stu_char(ii,mm);
        end
    end
    
% calculate the fraction for a given sample and then average over all samples
    CF.Wel_AA(1,pp,1) = mean(sum(wel(:,:,1))./ sum(MC.Stu_char));
    CF.Wel_AA(2,pp,1) = mean(sum(wel(:,:,2))./ sum(MC.Stu_char));
    CF.Wel_AA(3,pp,1) = mean(sum(wel(:,:,3))./ sum(MC.Stu_char));
    CF.Wel_AA(4,pp,1) = mean(sum(wel(:,:,4))./ sum(MC.Stu_char));

    CF.Wel_AA_sd(1,pp,1) = std(sum(wel(:,:,1))./ sum(MC.Stu_char));
    CF.Wel_AA_sd(2,pp,1) = std(sum(wel(:,:,2))./ sum(MC.Stu_char));
    CF.Wel_AA_sd(3,pp,1) = std(sum(wel(:,:,3))./ sum(MC.Stu_char));    
    CF.Wel_AA_sd(4,pp,1) = std(sum(wel(:,:,4))./ sum(MC.Stu_char));    
    
    % better off
    wel = zeros(I,M,4);
    for mm = 1:M
        for ii = 1:I
            wel(ii,mm,1) = (CFROL.U_ij(ii,squeeze(Mis.Assigned(mm,ii,1)),mm,pp) < CFROL.U_ij(ii,squeeze(CFROL.Assigned(mm,ii,pp)),mm,pp)) .* MC.Stu_char(ii,mm);
            wel(ii,mm,2) = (CFWTT.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) < CFWTT.U_ij(ii,squeeze(CFWTT.Assigned(mm,ii,pp)),mm,pp)) .* MC.Stu_char(ii,mm);
            wel(ii,mm,3) = (CFST.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) < CFST.U_ij(ii,squeeze(CFST.Assigned(mm,ii,pp)),mm,pp)) .* MC.Stu_char(ii,mm);
            wel(ii,mm,4) = (CFMis.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) < CFMis.U_ij(ii,squeeze(CFMis.Assigned(mm,ii,pp)),mm,pp)) .* MC.Stu_char(ii,mm);
        end
    end
    CF.Wel_AA(1,pp,2) = sum(sum(wel(:,:,1))) / sum(sum(MC.Stu_char));
    CF.Wel_AA(2,pp,2) = sum(sum(wel(:,:,2))) / sum(sum(MC.Stu_char));
    CF.Wel_AA(3,pp,2) = sum(sum(wel(:,:,3))) / sum(sum(MC.Stu_char));
    CF.Wel_AA(4,pp,2) = sum(sum(wel(:,:,4))) / sum(sum(MC.Stu_char));
    
% calculate the fraction for a given sample and then average over all samples
    CF.Wel_AA(1,pp,2) = mean(sum(wel(:,:,1))./ sum(MC.Stu_char));
    CF.Wel_AA(2,pp,2) = mean(sum(wel(:,:,2))./ sum(MC.Stu_char));
    CF.Wel_AA(3,pp,2) = mean(sum(wel(:,:,3))./ sum(MC.Stu_char));
    CF.Wel_AA(4,pp,2) = mean(sum(wel(:,:,4))./ sum(MC.Stu_char));

    CF.Wel_AA_sd(1,pp,2) = std(sum(wel(:,:,1))./ sum(MC.Stu_char));
    CF.Wel_AA_sd(2,pp,2) = std(sum(wel(:,:,2))./ sum(MC.Stu_char));
    CF.Wel_AA_sd(3,pp,2) = std(sum(wel(:,:,3))./ sum(MC.Stu_char));    
    CF.Wel_AA_sd(4,pp,2) = std(sum(wel(:,:,4))./ sum(MC.Stu_char));    

    % no change
    wel = zeros(I,M,4);
    for mm = 1:M
        for ii = 1:I
            wel(ii,mm,1) = (CFROL.U_ij(ii,squeeze(Mis.Assigned(mm,ii,1)),mm,pp) == CFROL.U_ij(ii,squeeze(CFROL.Assigned(mm,ii,pp)),mm,pp)) .* MC.Stu_char(ii,mm);
            wel(ii,mm,2) = (CFWTT.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) == CFWTT.U_ij(ii,squeeze(CFWTT.Assigned(mm,ii,pp)),mm,pp)) .* MC.Stu_char(ii,mm);
            wel(ii,mm,3) = (CFST.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) == CFST.U_ij(ii,squeeze(CFST.Assigned(mm,ii,pp)),mm,pp)) .* MC.Stu_char(ii,mm);
            wel(ii,mm,4) = (CFMis.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) == CFMis.U_ij(ii,squeeze(CFMis.Assigned(mm,ii,pp)),mm,pp)) .* MC.Stu_char(ii,mm);
        end
    end
    CF.Wel_AA(1,pp,3) = sum(sum(wel(:,:,1))) / sum(sum(MC.Stu_char));
    CF.Wel_AA(2,pp,3) = sum(sum(wel(:,:,2))) / sum(sum(MC.Stu_char));
    CF.Wel_AA(3,pp,3) = sum(sum(wel(:,:,3))) / sum(sum(MC.Stu_char));
    CF.Wel_AA(4,pp,3) = sum(sum(wel(:,:,4))) / sum(sum(MC.Stu_char));
    
% calculate the fraction for a given sample and then average over all samples
    CF.Wel_AA(1,pp,3) = mean(sum(wel(:,:,1))./ sum(MC.Stu_char));
    CF.Wel_AA(2,pp,3) = mean(sum(wel(:,:,2))./ sum(MC.Stu_char));
    CF.Wel_AA(3,pp,3) = mean(sum(wel(:,:,3))./ sum(MC.Stu_char));
    CF.Wel_AA(4,pp,3) = mean(sum(wel(:,:,4))./ sum(MC.Stu_char));

    CF.Wel_AA_sd(1,pp,3) = std(sum(wel(:,:,1))./ sum(MC.Stu_char));
    CF.Wel_AA_sd(2,pp,3) = std(sum(wel(:,:,2))./ sum(MC.Stu_char));
    CF.Wel_AA_sd(3,pp,3) = std(sum(wel(:,:,3))./ sum(MC.Stu_char)); 
    CF.Wel_AA_sd(4,pp,3) = std(sum(wel(:,:,4))./ sum(MC.Stu_char)); 
    
end

CF.Wel_nAA = NaN(4,P,3); % for the non-minority group
CF.Wel_nAA_sd = NaN(4,P,3); % standard deviation
for pp = 1:P
    % worse off
    wel = zeros(I,M,4);
    for mm = 1:M
        for ii = 1:I
            wel(ii,mm,1) = (CFROL.U_ij(ii,squeeze(Mis.Assigned(mm,ii,1)),mm,pp) > CFROL.U_ij(ii,squeeze(CFROL.Assigned(mm,ii,pp)),mm,pp)) .* (1-MC.Stu_char(ii,mm));
            wel(ii,mm,2) = (CFWTT.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) > CFWTT.U_ij(ii,squeeze(CFWTT.Assigned(mm,ii,pp)),mm,pp)) .* (1-MC.Stu_char(ii,mm));
            wel(ii,mm,3) = (CFST.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) > CFST.U_ij(ii,squeeze(CFST.Assigned(mm,ii,pp)),mm,pp)) .* (1-MC.Stu_char(ii,mm));
            wel(ii,mm,4) = (CFMis.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) > CFMis.U_ij(ii,squeeze(CFMis.Assigned(mm,ii,pp)),mm,pp)) .* (1-MC.Stu_char(ii,mm));
        end
    end
% calculate the fraction for a given sample and then average over all samples
    CF.Wel_nAA(1,pp,1) = mean(sum(wel(:,:,1))./ sum(1-MC.Stu_char));
    CF.Wel_nAA(2,pp,1) = mean(sum(wel(:,:,2))./ sum(1-MC.Stu_char));
    CF.Wel_nAA(3,pp,1) = mean(sum(wel(:,:,3))./ sum(1-MC.Stu_char));
    CF.Wel_nAA(4,pp,1) = mean(sum(wel(:,:,4))./ sum(1-MC.Stu_char));

    CF.Wel_nAA_sd(1,pp,1) = std(sum(wel(:,:,1))./ sum(1-MC.Stu_char));
    CF.Wel_nAA_sd(2,pp,1) = std(sum(wel(:,:,2))./ sum(1-MC.Stu_char));
    CF.Wel_nAA_sd(3,pp,1) = std(sum(wel(:,:,3))./ sum(1-MC.Stu_char));
    CF.Wel_nAA_sd(4,pp,1) = std(sum(wel(:,:,4))./ sum(1-MC.Stu_char));
    
    % better off
    wel = zeros(I,M,4);
    for mm = 1:M
        for ii = 1:I
            wel(ii,mm,1) = (CFROL.U_ij(ii,squeeze(Mis.Assigned(mm,ii,1)),mm,pp) < CFROL.U_ij(ii,squeeze(CFROL.Assigned(mm,ii,pp)),mm,pp)) .* (1-MC.Stu_char(ii,mm));
            wel(ii,mm,2) = (CFWTT.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) < CFWTT.U_ij(ii,squeeze(CFWTT.Assigned(mm,ii,pp)),mm,pp)) .* (1-MC.Stu_char(ii,mm));
            wel(ii,mm,3) = (CFST.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) < CFST.U_ij(ii,squeeze(CFST.Assigned(mm,ii,pp)),mm,pp)) .* (1-MC.Stu_char(ii,mm));
            wel(ii,mm,4) = (CFMis.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) < CFMis.U_ij(ii,squeeze(CFMis.Assigned(mm,ii,pp)),mm,pp)) .* (1-MC.Stu_char(ii,mm));
        end
    end
% calculate the fraction for a given sample and then average over all samples
    CF.Wel_nAA(1,pp,2) = mean(sum(wel(:,:,1))./ sum(1-MC.Stu_char));
    CF.Wel_nAA(2,pp,2) = mean(sum(wel(:,:,2))./ sum(1-MC.Stu_char));
    CF.Wel_nAA(3,pp,2) = mean(sum(wel(:,:,3))./ sum(1-MC.Stu_char));
    CF.Wel_nAA(4,pp,2) = mean(sum(wel(:,:,4))./ sum(1-MC.Stu_char));

    CF.Wel_nAA_sd(1,pp,2) = std(sum(wel(:,:,1))./ sum(1-MC.Stu_char));
    CF.Wel_nAA_sd(2,pp,2) = std(sum(wel(:,:,2))./ sum(1-MC.Stu_char));
    CF.Wel_nAA_sd(3,pp,2) = std(sum(wel(:,:,3))./ sum(1-MC.Stu_char));
    CF.Wel_nAA_sd(4,pp,2) = std(sum(wel(:,:,4))./ sum(1-MC.Stu_char));

    % no change 
    wel = zeros(I,M,4);
    for mm = 1:M
        for ii = 1:I
            wel(ii,mm,1) = (CFROL.U_ij(ii,squeeze(Mis.Assigned(mm,ii,1)),mm,pp) == CFROL.U_ij(ii,squeeze(CFROL.Assigned(mm,ii,pp)),mm,pp)) .* (1-MC.Stu_char(ii,mm));
            wel(ii,mm,2) = (CFWTT.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) == CFWTT.U_ij(ii,squeeze(CFWTT.Assigned(mm,ii,pp)),mm,pp)) .* (1-MC.Stu_char(ii,mm));
            wel(ii,mm,3) = (CFST.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) == CFST.U_ij(ii,squeeze(CFST.Assigned(mm,ii,pp)),mm,pp)) .* (1-MC.Stu_char(ii,mm));
            wel(ii,mm,4) = (CFMis.U_ij(ii,squeeze(Mis.Assigned(mm,ii,pp)),mm,pp) == CFMis.U_ij(ii,squeeze(CFMis.Assigned(mm,ii,pp)),mm,pp)) .* (1-MC.Stu_char(ii,mm));
        end
    end
% calculate the fraction for a given sample and then average over all samples
    CF.Wel_nAA(1,pp,3) = mean(sum(wel(:,:,1))./ sum(1-MC.Stu_char));
    CF.Wel_nAA(2,pp,3) = mean(sum(wel(:,:,2))./ sum(1-MC.Stu_char));
    CF.Wel_nAA(3,pp,3) = mean(sum(wel(:,:,3))./ sum(1-MC.Stu_char));
    CF.Wel_nAA(4,pp,3) = mean(sum(wel(:,:,4))./ sum(1-MC.Stu_char));

    CF.Wel_nAA_sd(1,pp,3) = std(sum(wel(:,:,1))./ sum(1-MC.Stu_char));
    CF.Wel_nAA_sd(2,pp,3) = std(sum(wel(:,:,2))./ sum(1-MC.Stu_char));
    CF.Wel_nAA_sd(3,pp,3) = std(sum(wel(:,:,3))./ sum(1-MC.Stu_char));    
    CF.Wel_nAA_sd(4,pp,3) = std(sum(wel(:,:,4))./ sum(1-MC.Stu_char));    
end

clear wel Mis.Assigned CFROL.Assigned CRTT.Assigned CRST.Assigned U_ij

CF.diff_wel = CF.Wel(:,:,2) - CF.Wel(:,:,1); % 4xP
CF.diff_wel_AA = CF.Wel_AA(:,:,2) - CF.Wel_AA(:,:,1);
CF.diff_wel_nAA = CF.Wel_nAA(:,:,2) - CF.Wel_nAA(:,:,1);

CF.Wel_AA = CF.Wel_AA* 100;
CF.Wel_nAA = CF.Wel_nAA* 100;
CF.Wel_AA_sd = CF.Wel_AA_sd* 100;
CF.Wel_nAA_sd = CF.Wel_nAA_sd* 100;

% --------------- %
% Table: Welfare  %
% --------------- %

tab_wel_AA = table('Size',[5*P-1,8], ...
    'VariableTypes',{'double','double','double','double','double','double','double','double'},...
    'VariableNames',{'Worse off (mean)','Worse off (s.d.)','|','Better off (mean)','Better off (s.d.)',...
    '"','Indifferent (mean)', 'Indifferent (s.d.)'},...
    'RowNames',{'DGP TT: Submitted ROLs','DGP TT: WTT','DGP TT: Stability','DGP TT: Actual behavior (the truth)',...
    '-','DGP PIM: Submitted ROLs','DGP PIM: WTT','DGP PIM: Stability','DGP PIM: Actual behavior (the truth)',...
    '~','DGP PRM: Submitted ROLs','DGP PRM: WTT','DGP PRM: Stability','DGP PRM: Actual behavior (the truth)'});

% tab_wel_AA = NaN(5*P-1,6);
% AA
for pp = 1:P
     % worse off: mean
    tab_wel_AA((pp-1)*5+1:(pp-1)*5+4,1) = table(CF.Wel_AA(1:4,pp,1)); 

     % worse off: sd
    tab_wel_AA((pp-1)*5+1:(pp-1)*5+4,2) = table(CF.Wel_AA_sd(1:4,pp,1)); 
    
     % better off: mean
    tab_wel_AA((pp-1)*5+1:(pp-1)*5+4,4) = table(CF.Wel_AA(1:4,pp,2));
    
     % better off: sd
    tab_wel_AA((pp-1)*5+1:(pp-1)*5+4,5) = table(CF.Wel_AA_sd(1,pp,2)); 

     % no change: mean
    tab_wel_AA((pp-1)*5+1:(pp-1)*5+4,7) = table(CF.Wel_AA(1,pp,3)); 

     % no change: mean
    tab_wel_AA((pp-1)*5+1:(pp-1)*5+4,8) = table(CF.Wel_AA_sd(1,pp,3));           
end

% non-AA

tab_wel_nAA = table('Size',[5*P-1,8], ...
    'VariableTypes',{'double','double','double','double','double','double','double','double'},...
    'VariableNames',{'Worse off (mean)','Worse off (s.d.)','|','Better off (mean)','Better off (s.d.)',...
    '"','Indifferent (mean)', 'Indifferent (s.d.)'},...
    'RowNames',{'DGP TT: Submitted ROLs','DGP TT: WTT','DGP TT: Stability','DGP TT: Actual behavior (the truth)',...
    '-','DGP PIM: Submitted ROLs','DGP PIM: WTT','DGP PIM: Stability','DGP PIM: Actual behavior (the truth)',...
    '~','DGP PRM: Submitted ROLs','DGP PRM: WTT','DGP PRM: Stability','DGP PRM: Actual behavior (the truth)'});

%     tab_wel_nAA = NaN(5*P,6);
for pp = 1:P
     % worse off: mean
    tab_wel_nAA((pp-1)*5+1:(pp-1)*5+4,1) = table(CF.Wel_nAA(1:4,pp,1)); 

     % worse off: sd
    tab_wel_nAA((pp-1)*5+1:(pp-1)*5+4,2) = table(CF.Wel_nAA_sd(1:4,pp,1)); 
    
     % better off: mean
    tab_wel_nAA((pp-1)*5+1:(pp-1)*5+4,4) = table(CF.Wel_nAA(1:4,pp,2));
    
     % better off: sd
    tab_wel_nAA((pp-1)*5+1:(pp-1)*5+4,5) = table(CF.Wel_nAA_sd(1,pp,2)); 

     % no change: mean
    tab_wel_nAA((pp-1)*5+1:(pp-1)*5+4,7) = table(CF.Wel_nAA(1,pp,3)); 

     % no change: mean
    tab_wel_nAA((pp-1)*5+1:(pp-1)*5+4,8) = table(CF.Wel_nAA_sd(1,pp,3));           
end

cd(folder_tables)
writetable(tab_wel_AA,'tab_E2_wel_AA.csv','WriteRowNames',true)
writetable(tab_wel_nAA,'tab_E2_wel_nAA.csv','WriteRowNames',true)

tab_T = table('Size',[2,2],'VariableTypes',{'double','double'},'VariableNames',{'mean','std'},'RowNames',{'T1','T0'});

tab_T(1,1) = table(mean(sum(MC.Stu_char)));
tab_T(1,2) = table(std(sum(MC.Stu_char)));
tab_T(2,1) = table(mean(sum(1-MC.Stu_char)));
tab_T(2,2) = table(std(sum(1-MC.Stu_char)));

% figures

cd(folder_figures)

% ------ %
%   AA   %
% ------ %

figure()
% area([2;3;4],[1*ones(3,1)],'FaceColor',[0.9  0.9  0.9],'LineStyle','none');
hold on
x = 1:P;
p1 = plot(x,CF.diff_wel_AA(1,:),'-.g','LineWidth',2);
hold on
p2 = plot(x,CF.diff_wel_AA(2,:),'-r','LineWidth',2);
hold on
p3 = plot(x,CF.diff_wel_AA(3,:),':b','LineWidth',2);
hold on
p4 = plot(x,CF.diff_wel_AA(4,:),'--m','LineWidth',2);
hold on
%line([1 8], [CF.diff_wel_AA(1,1) CF.diff_wel_AA(1,1)],'Color','black','LineWidth',0.5)
legend([p1,p2,p3,p4],{'Estimates w/ Submitted ROLs','WTT Estimates', 'Stability Estimates','Actual Behavior (the truth)'},'Box','off','Location','southwest')

% x0 = [2/8 1.1/8];
% y0 = [0.85 0.8];
% annotation('textarrow',x0,y0,'String','true value (based on true preferences) ')

%
% x1 = 2.25;
% y1 = 0.5;
% txt1 = 'Payoff Irrelevant';
% text(x1,y1,txt1)
%
% x2 = 5.75;
% y2 = 0.5;
% txt2 = 'Payoff Relevant';
% text(x2,y2,txt2)

%title(figname)
ax = gca; % current axes
ax.FontSize = 14;
ax.TickDir = 'out';
ax.TickLength = [0.01 0.01];
ax.YLim = [0.7 1];
% ax.XLim = [1 4];
% ax.XTick = ([1 2 3 4]);
% ax.XTickLabel = ({'TT','IRR1','IRR2','REL'});
ax.XLim = [1 3];
ax.XTick = ([1 2 3]);
ax.XTickLabel = ({'TT','PIM','PRM'});
xlabel('Data Generating Processes');
ylabel('( Fraction of winners ) - ( Fraction of losers )');
print('fig_4_CF_wel_AA','-djpeg')
print('fig_4_CF_wel_AA','-depsc2')

% ------ %
% non-AA %
% ------ %

figure()
% area([2;3;4],-ones(3,1),'FaceColor',[0.9  0.9  0.9],'LineStyle','none');
hold on
x = 1:P;
p1 = plot(x,CF.diff_wel_nAA(1,:),'-.g','LineWidth',2);
hold on
p2 = plot(x,CF.diff_wel_nAA(2,:),'-r','LineWidth',2);
hold on
p3 = plot(x,CF.diff_wel_nAA(3,:),':b','LineWidth',2);
% xticks([1 2 3 4 5 6 7 8])
% xticklabels({'TT','payoff','irrelevant',' ','skips','payoff','relevant','mistakes'})
hold on
p4 = plot(x,CF.diff_wel_nAA(4,:),'--m','LineWidth',2);
hold on
legend([p1,p2,p3,p4],{'Estimates w/ Submitted ROLs','WTT Estimates', 'Stability Estimates','Actual Behavior (the truth)'},'Box','off','Location','northwest')

% line([1 8], [CF.diff_wel_nAA(1,1) CF.diff_wel_nAA(1,1)],'Color','black','LineWidth',0.5)

% x1 = 1.5;
% y1 = -0.65;
% txt1 = '\uparrow the true value (based on true preferences)';
% text(x1,y1,txt1)

% x1 = 2.25;
% y1 = -0.2;
% txt1 = 'Payoff Irrelevant';
% text(x1,y1,txt1)
%
% x2 = 5.75;
% y2 = -0.2;
% txt2 = 'Payoff Relevant';
% text(x2,y2,txt2)

% x0 = [2/8 1.1/8];
% y0 = [0.375 0.45];
% annotation('textarrow',x0,y0,'String','true value (based on true preferences) ')

ax = gca; % current axes
ax.FontSize = 14;
ax.TickDir = 'out';
ax.TickLength = [0.01 0.01];
ax.YLim = [-0.7 -0.4];
% ax.XLim = [1 4];
% ax.XTick = ([1 2 3 4]);
% ax.XTickLabel = ({'TT','IRR1','IRR2','REL'});
ax.XLim = [1 3];
ax.XTick = ([1 2 3]);
ax.XTickLabel = ({'TT','PIM','PRM'});
xlabel('Data Generating Processes');
ylabel('( Fraction of winners ) - ( Fraction of losers )');
print('fig_E2_CF_wel_nAA','-djpeg')
print('fig_E2_CF_wel_nAA','-depsc2')

close all

toc